Equilibration times in numerical simulation of structural glasses: Comparing parallel 
tempering and conventional molecular dynamics 
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Generation of equilibrium configurations is the major obstacle for numerical investigation of the 
slow dynamics in supercooled liquid states. The parallel tempering (PT) technique, originally pro- 
posed for the numerical equilibration of discrete spin-glass model configurations, has recently been 
applied in the study of supercooled structural glasses. We present an investigation of the ability 
of parallel tempering to properly sample the liquid configuration space at different temperatures, 
by mapping the PT dynamics into the dynamics of the closest local potential energy minima (in- 
herent structures). Comparing the PT equilibration process with the standard molecular dynamics 
equilibration process we find that the PT does not increase the speed of equilibration of the (slow) 
configurational degrees of freedom. 

PACS numbers: 61.20.Ja, 02.70.Ns, 02.70.Uu, 61.43.Fs, 64.70.Pf 



As a liquid is cooled below its melting temperature T m 
(supercooled liquid) the structural time r increases con- 
siderably. In a small temperature interval, t changes by 
more than 13 order of magnitude. When r reaches val- 
ues bigger than 100s the liquid behaves as an amorphous 
solid, i.e. a glass. 

In recent years, a considerable interest has been de- 
voted to the study of the supercooled state of matter, 
both theoretically H, p[ [| , experimentally M, [|,J| and nu- 
merically (7[ H [lQP Both thermodynamics |gf and dy- 
namic [yj theories have been proposed to explain the rich 
phenomenology of glassy systems. Molecular dynamics 
(MD) simulations have proved to be a powerful tool for 
studying simple models for liquids in supercooled states 
(for a review, see Jll[]). Simulation stretching in the ns 
time window have offered the possibility of a detailed 
comparison between theoretical predictions and "exact" 
numerical results. So far, such comparisons have been 
limited to weakly supercooled states, i.e. to the temper- 
ature region where characteristic times are at most of 
the order of 10 ns. In this region, mode coupling the- 
ory (MCT) has shown its ability in correctly predicting 
the numerical results |j| |lj, [l3| even for network forming 
liquids @ ||. 

The analysis of numerical data has been also very fruit- 
ful in the study of the potential energy surface (PES) - 
the so-called energy landscape — of several models [||. 
These studies have provided evidence that in equilibrium 
the average basin depth ejs(T) is a decreasing function of 
T |l6| . The number of explored local PES minima, com- 
monly named inherent structures (IS) — the exponential 
of the configurational entropy in the inherent structure 
formalism |l7], |l8[ |l^] — decreases also on cooling. Nu- 
merical studies on aging liquids |2(], |2lJ have shown that 
the equilibration process is related to the slow search for 
deeper and deeper basins on the potential energy surface. 
In the PES framework at least two different factors con- 



trol the equilibration time scale: (i) the timescale for es- 
cape from a selected basin (a timescale depending on the 
kinetic energy) and (ii) the timescale for finding deeper 
basins (a time scale depending on the number of accessi- 
ble basins). Which of the two factors is the leading one 
is still an open question. 



Presently, the interesting region where dynamics slow 
down beyond the ns time scale can not be studied numer- 
ically since the generation of equilibrated configurations 
requires prohibitive computational times. The possibil- 
ity of disposing of equilibrium configuration could open 
the possibility of studying, if not the entire structural 
relaxation process, at least the initial part of it, where 
several interesting phenomena related to the dynamics 
in disordered structures are taking place [E2l E3l E4l E5|. 



Several algorithms have been developed to improve the 
equilibration times in numerical simulations of glassy sys- 
tems |^6|, I!?! . A study by Kob and Yamamoto suggests 
that the parallel tempering (PT) may become an impor- 
tant tool to provide independent equilibrium configura- 
tions for structural glasses. The PT technique [g7| was 
developed for dealing with the slow dynamics of disor- 
dered spin systems. The PT algorithms simultaneously 
simulates a set of M identical non-interacting replicas of 
the system, each of them at a different T. Pairs of replica 
swap their temperatures according to a Monte Carlo pro- 
cedure. The basic idea is that each replica performs a 
random walk among the M different T. Hence, when 
the replica explores the high T states, the probability to 
escape from its basin is enhanced. 



Each of the M replicas, composed by N atoms, is de- 
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scribed by an Hamiltonian 
N 1 

H m (q m ,p m ) = ^ + A m (t)E(q m ) + 

f — ' Am 

where E(q m ) is the potential energy of the system. A m (t) 
is a scaling parameter for the potential energy, which 
effectively sets the temperature T of the m-th replica to 
the value Xb/A m (i), where To indicates the lowest studied 
temperature. Consequently the values A m (0), for m — 
0....M — 1 set the M different temperatures of the M 
replicas |@ at time 0. The degree of freedom s m (last two 
terms in Eq.|l|) are relative to the Nose thermal bath [^8| . 
The thermostat constrains the average kinetic energy of 
each replica to the value 3/2A^fcsTo. 
The whole Hamiltonian is then: 

M 

H=J2 H ™- (2) 

As discussed in detail in Ref. |27j] and Q , the choice of 
A m (0) must guarantee a significant overlap in the energy 
distributions of different replicas, a requirement which 
oblige to keep M proportionally to the system size. 

In this article we focus on the time requested to find the 
low inherent structure configurations visited in equilib- 
rium. More specifically, to evaluate if the PT technique 
is a viable candidate to equilibrate structural glasses, we 
compare the PT and the conventional MD dynamics by 
computing the inherent structure energy as a function of 
the simulation time. Since e/s is a much more sensitive 
indicator of equilibrium than the total potential energy, 
we can put the PT technique under stringent test. 

I. MODELS AND DETAILS OF SIMULATION 

The system we investigated is the monoatomic 
Lennard-Jones (LJ) model modified by adding a many- 
body anti-crystalline potential designed to inhibit crys- 
tallization The e and a parameter of the LJ poten- 
tial are chosen as unit of energy and length respectively. 
The LJ potential is truncated and shift at 2.5. The poten- 
tial energy, which includes the anti-crystalline potential 
is 

E(q m ) = V LJ {q m ) + 

ia^0(S m (fc)-S o )[S m (fc)-S o ] 2 . (3) 

% 

where Vlj is the LJ part of the potential and the sum 
is over all k such that k max — Ak < \\k\\ < k max + Ak. 
The other terms in Eq.Q set in only when any of the 
wavevector around the structure factor peak increases 



beyond the value So and act by damping the unwanted 
crystalline like density fluctuation. We chose a number 
density p — 1 , and So — 10, k max — 7.12, a = 0.83 and 
Afc = 0.34 for the anti-crystalline parameters as proposed 
in Ref. fji] . The integration timestep is 0.0025, in time 
units of yjmo 2 /e. The dynamics for this model has been 
previously studied |52). It has been shown that a fast 
decrease of the structural times takes place below T = 1 . 
The T dependence of r follows a power law in T — T x , 
with T x w 0.475. T x has been identified with the ideal 
MCT for this model, an hypothesis supported also by an 
analysis of the T dependence of the diffusive directions 

The PT algorithm is identical to that encoded in 
Ref. (3^] and we refer to that paper for details on the tech- 
nique. The algorithm we implement uses M — 14 iden- 
tical non-interacting replicas each composed of TV = 256 
particles. The 14 temperatures are chosen to span an 
range from T — 1.05 down to 0.485, in particular the 
temperatures we used are the following: 0.485, 0.518, 
0.534, 0.562, 0.597, 0.646, 0.694, 0.745, 0.80, 0.85, 0.90, 
0.95, 1.0, 1.05. 

The Hamiltonian of one replica is that of Eq. [[[) . All 
replicas evolve according to the standard Nose' constant 
temperature MD simulation. Every 1000 steps an at- 
tempt to exchange the scaling parameter of all pair of 
replicas with adjacent temperatures (swap of the A val- 
ues) is performed using the following criterion: an ex- 
change is accepted in Metropolis fashion, i.e. the accep- 
tance ratio is: 

A m ,„ < 

Wm ' n ~ \ exp(-A m ,„) , A m ,„ > 

where A TOi „ = (3 (A n - A m )[E(q m ) - E(q n )]. The events 
with i = 0, 2, 4, . . . or i = 1, 3, 5, . . . are repeated alter- 
natively every 1000 integration steps. 

The outcome of such calculation are, in principle, equi- 
librium configurations in the canonical ensemble at the 
M different temperatures. 

To estimate the time requested to the PT algorithm 
to equilibrated all replicas we start our PT algorithm 
with M replicas extracted from a previously generated 
ensemble of equilibrium configurations at T = 1.05. At 
this T, the structural relaxation time is of the order of 
1000 steps and hence generation of equilibrium configura- 
tions with conventional MD does not pose any problem. 
By starting with this ensemble of configurations, the PT 
equilibration time is by construction for the highest 
temperature. We performed 20 of such independent PT 
runs to improve the statistic. Each of such run lasts two 
million time step. Hence, in the PT part of the work, 
the equation of motion have been integrated more than 
40 million times. 

The same starting configurations (T = 1.05) are also 
used as initial configurations for conventional constant 
temperature MD simulations at several bath tempera- 
tures T, to compare the rate of equilibration of PT and 
conventional MD algorithms. For each temperature of 
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these MD simulations we performed 16 independent runs 
to improve statistics. 

Local minima configurations have been calculated via 
conjugate gradient minimization. The minimization pro- 
cess is considered completed when the potential energy 
change associated with one iteration is less than 10 -15 to 
ensure a great accuracy. 



II. RESULTS 

We focus here on the evolution in time of the inherent 
structure energy e/s, comparing the PT and MD proce- 
dures. Recent work [^0[ |2l], |33) has provided evidence 
that following a T change, the system response is charac- 
terized by two different time scales. A short time scale, 
related to the equilibration of the system within a well 
defined basin of the energy landscape and a slow time 
scale related to the search for basin of the "right" depth. 
The evolution of one-time quantities carries on informa- 
tion on these two time scales. After a sudden change of 
T, a very fast decreases of the vibrational energy — cor- 
responding to the fast equilibration to the new bath T of 
the intra-basin vibrational motions — is observed. This 
fast change is followed by a much slower decrease which 
corresponds to the slow decrease of the basin's depth. 
The absolute change in e/s during the aging process is 
significantly smaller than the change in the total poten- 
tial energy and hence it requires a careful analysis to be 
detected. A way which has been proven fruitful to sep- 
arate large fast component and small slow component is 
to monitor directly the evolution of e/s- Building on the 
expertise developed in recent years, we adopt this indi- 
cator as effective tool to monitor the equilibration of the 
system in configuration space. 

We note on passing that, since e/s is a small com- 
ponent of the potential energy — being the intra-basin 
vibrational part dominating — a nice scaling of the total 
potential energy distribution may not guarantee perfect 
equilibration. 

Figure |l| compares the time evolution of the inherent 
structure energy in the conventional MD (Top) and in the 
PT (Bottom) runs. In both cases, by construction, the 
initial e/s coincides with the equilibrium value of e/s at 
T = 1.05. As time goes on, each replica starts to explore 
larger and larger parts of the configuration space select- 
ing configuration with lower and lower values of e/s- The 
equilibration process lasts until the equilibrium value of 
e/s is reached. The same picture applies to the conven- 
tional MD case, where the simulation indeed reproduces 
the aging process following a T jump from T = 1.05 to 
the new bath temperature. 

Fig. U compares the time evolution of e/s for PT and 
MD simulations at three different temperatures. In all 
cases, we find clear indications that the equilibration of 
the slow degrees of freedom does not depend on the pro- 
cedure adopted. 
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FIG. 1: Inherent structure energies as a function of time for 
MD (Top) and PT (Bottom). 
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FIG. 2: Inherent structure energies as a function of time: 
comparison between MD and PT . 
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III. DISCUSSION 

The data shown in Fig.^J very clearly show that for 
the Lennard Jones case investigated in this article no 
improvement in equilibration rates is achieved by imple- 
menting the PT algorithm. 

For the related model of binary mixture Lennard Jones 
p| , the number of basins with depth ejs, in the range 
of eis values characteristic of the PES region explored 
above T x is well represented by a gaussian distribution. 
The total number of basins has been shown to scale with 
the size N of the system as e aN with a « 0.8 pl|. For 
our 256 atom case, this correspond to about 10 8 ^basins. 
Under such complicated potential energy landscape con- 
ditions, an unbiased search for the location of the deepest 
basin would require an order of 10 87 attempts! In this 
respect, it is possible that the rate of equilibration at 
low temperatures is significantly controlled by a simple 
entropic effect. This could explain why the possibility 
of overcoming barriers with higher probability offered by 
PT does not favor a faster equilibration process. This pic- 
ture is also consistent with the fact that in the T-region 
explored (which is still above T x ) saddle dominated dy- 
namics is dominant. Recent instantaneous normal modes 
analysis |38| has indeed provided evidence that above T x 
the system explore mostly regions of the potential energy 
landscape which are characterized by a large number of 
negative curvature directions. No activated processes are 
requested in this condition to change local basin. 

It would be interesting to find out if the PT algorithm 
may be valuable in studying strong liquids, for which less 
relevant changes in the PES are taking place on cooling 
as compared to fragile liquids Q and for which activated 
pro cesses are dominant at low T . Preliminary indications 
P5| seems to suggest that this may be the case. It would 
also be important to correlate the efficiency of PT with 
the structure of configuration space and connectivity be- 



tween distinct potential energy surface basins. A possible 
line of research could be to compare eis{t) for PT and 
MD in clusters with different disconnectivity graph |ol| 
types. 

The fact that in the explored T-region equilibration 
times are not improved in the PT case is consistent with 
recent finding in the field of protein simulations. This 
analogy is not at all surprising due to the significant sim- 
ilarities in the problem of glass formation and protein 
folding p7| ]. Again, the similarity points to the structure 
of the PES as key element in the control of the charac- 
teristic times. 

To conclude, we like to call the reader attention on 
the fact that to study a fixed T-range, the PT tech- 
nique requires a number of replicas which increases lin- 
early with the system size, to guarantee proper overlap 
in the potential energy distributions of adjacent replicas 
and hence a significant replica exchange rate. Moreover, 
all replicas have to be simulated for the same total time 
interval. This time is fixed by the lowest temperature, 
which is characterized by a relaxation time which may 
well be several order of magnitude smaller than the one 
of the top temperature. Both effects concur in making 
PT a not very convenient algorithm for simulating struc- 
tural glasses as compared to conventional MD. Indeed, in 
MD, the total simulation time at each temperature can 
be chosen to scale with the structural relaxation time. 
Of course, bookkeeping facilities are enhanced in the PT 
case. 
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